An integral with a square root of a polynomial in its denominator

$I\left(x\right)=\int \frac{dx}{\sqrt{P\left(x\right)}}$

can be expressed as an inverse sine function if the polynomial is quadratic. For cubic or quartic polynomials it is one of the elliptic integrals, which are inverse elliptic functions. If the degree of the polynomial is greater than four it is called a hyperelliptic integral and is not generally expressible in terms of known special functions.

The word “polynomial” as used by mathematicians means integer exponents only. This ensures that polynomials obey addition and multiplication as abstracted in rings, making available the mathematical apparatus of abstract algebra. It also ensures that hyperelliptic integrals are single-valued functions in the complex plane.

Integration of the classical energy equation for power potentials with arbitrary real exponent can lead to an integral of the form

$I(x,k)=\int \frac{dx}{\sqrt{1-x^k}}$

Since the exponent k is allowed to be nonintegral this integral is not a standard hyperelliptic integral. It is multiple-valued as a complex function and does not have any nice symmetry properties that allow application of known techniques for inversion of hyperelliptic integrals. One nevertheless would like an expression for its inverse on the principal Riemann sheet.

The generalized hyperelliptic integral can always be evaluated with a Gauss hypergeometric function

$I(x,k)=x{}_{2}F_1[\frac{1}{2},\frac{1}{k};\frac{1}{k}+1;x^k]$

For k = 2 this hypergeometric function is equivalent to an inverse sine as expected, but will be left in this form to express the integral for all positive real values of the exponent. Negative real exponent values require changes to avoid a singular argument and a related imaginary unit, and will be considered later.

For application to classical motions, the generalized hyperelliptic integral is equated to a temporal variable

$I\left[x\right(t),k]=t$

and the explicit temporal dependence of the motion is found by inverting this equation. The inversion can be easily visualized with a parametric plot whose abscissa is the generalized hyperelliptic integral and whose ordinate is the independent variable:

This visual inversion represents the first quarter of an oscillatory motion whose period is determined by the singularity of the integrand at x = 1 . The value of the integral between zero and and that point is

$\frac{T}{4}=\frac{\Gamma (\frac{1}{k}+1)\Gamma \left(\frac{1}{2}\right)}{\Gamma (\frac{1}{k}+\frac{1}{2})\Gamma \left(1\right)}=\frac{\sqrt{\pi }\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}$

using an expression for the value of the hypergeometric function at unity. The full period can be written equivalently

$T=4\sqrt{\pi }\frac{\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}=\frac{2^{1+2/k}}{k}\frac{\Gamma (\frac{1}{k}{)}^2}{\Gamma \left(\frac{2}{k}\right)}$

using the Legendre duplication formula. For k = 2 the full period is 2π as expected.

One way to achieve inversion of the generalized hyperelliptic integral is with Lagrange inversion: for a function $t=f\left(x\right)$ , the inverse function expanded about a center c is

$x=f^{-1}\left(t\right)=c+\sum _{m=1}^{\infty }\frac{[t-f(c){]}^m}{m!}\underset{x\to c}{lim}\frac{d^{m-1}}{d{x}^{m-1}}[\frac{x-c}{f\left(x\right)-f\left(c\right)}{]}^m$

The limits in this expression of multiple derivatives of powers of the rational function can be evaluated using Mathematica or its equivalent for integral values of the exponent k, but Mathematica at least does not yet handle fractional calculus. Unfortunately the evaluation rapidly becomes too involved to yield a practical result for large integral exponents.

What one can get from evaluation of the limits is an indication of the general form of the power series that represents an inversion of the generalized hyperelliptic integral. Knowing the appropriate form of the expansion allows one to arrange series expansions in ways that Mathematica can handle readily, thereby automating the determination of expansion coefficients. Many of these evaluations take some time to run, so only final results will be presented.

For motion that begins at the origin there are two centers of expansion that are fruitful, the origin x = 0 and the turning point x = 1 . Motion that begins somewhere between the origin and the turning point can be expressed by shifting the temporal variable appropriately.

The limits of multiple derivatives of powers of the rational function for an expansion about the origin present no difficulties, and they are

The limits for an expansion about the turning point are not so simply determined. Using an exact limit of unity produces imaginary singularities for limits of odd-order multiple derivatives, in addition to values that are half of what is expected: presumably this occurs because the singularity in the integral falls on the radius of convergence of the series and Mathematica assigns only half of the value of complex residues determining the results. The expected values can be deduced by using a rational approximation to unity that indicates approximate limits, but the result is not clean as can be seen in the code. The way to get a clean answer is by adding a factor of two for the exact limit of unity and removing imaginary values explicitly. The correct values are

It is worth noting that simply graphing the functions resulting from the multiple derivatives will not indicate the correct limits for the expansion about the turning point. Limit is more careful than Plot in combining large quantities that ultimately cancel one another.

The results for limits indicate that the expansion about the origin is of the form

$x_{\mathrm{origin}}\left(t\right)=t+\sum _{m=1}^{\infty }(-1{)}^m{c}_m{t}^{mk+1}$

while the expansion about the turning point is of the form

$x_{\mathrm{turning}}\left(t\right)=1+\sum _{m=1}^{\infty }(-1{)}^m{d}_m[t-\frac{\sqrt{\pi }\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}{]}^{2m}$

Another way to justify the form of the expansion about the origin is to set $x_{\mathrm{origin}}\left(t\right)=tf\left(t\right)$ in the equation to be inverted and cancel a single power of the temporal variable,

$f\left(t\right){}_{2}F_1[\frac{1}{2},\frac{1}{k};\frac{1}{k}+1;t^k\left[f\right(t){]}^k]=1$

and it is immediately clear that $f\left(t\right)$ can only be a function of powers of $t^k$ . The form of the expansion about the turning point is a generalization of the behavior of sine about $\frac{\pi }{2}$ , which is equivalent to a cosine about the origin.

The expected procedure at this point is to expand the generalized hyperelliptic integral about either point in Mathematica, insert the expansion, expand and collect on equal powers and set equal to the temporal variable to determine coefficients recursively. Unfortunately Mathematica (along with humans) will balk at arbitrary powers of an infinite summation. It is much more convenient to work with expansions in product form,

$\begin{array}{l}{x}_{\mathrm{origin}}\left(t\right)=t\prod _{m=1}^{\infty }(1+a_m{t}^{mk})\\ x_{\mathrm{turning}}\left(t\right)=\prod _{m=1}^{\infty }(1+b_m{u}^{2m})u=t-\frac{\sqrt{\pi }\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}\end{array}$

because then the arbitrary powers of binomials can be expanded in general form as

$(1+y{)}^k=\sum _m\left(\genfrac{}{}{0ex}{}{k}{m}\right)y^m=\sum _m\frac{\Gamma (k+1)}{\Gamma (k-m+1)m!}{y}^m$

which is how Mathematica interprets binomial coefficients of nonintegral argument in any case.

For the expansion about the origin, the nonintegral power in the argument of the hypergeometric function will lead to singular expansion coefficients if the power series is developed with respect to x directly. Instead since a power is small whenever its base is small, the power series is developed with respect to the entire argument x^k and the arbitrary power substituted afterward. The justification for this process is that it leads to a sensible result that works.

Mathematica then makes relatively light work of integral powers of series to the desired level of accuracy, from which terms can be collected and groups of coefficients equated to zero. The result for an expansion depth of four terms beyond the initial linear is

$\begin{array}{l}{x}_{\mathrm{origin}}\left(t\right)=t[1-\frac{1}{2(k+1)}{t}^k][1+\frac{(k-1)}{8(k+1)(2k+1)}{t}^{2k}]\\ \times [1-\frac{(k-1)(k^2-4k-2)}{16(k+1{)}^2(2k+1\left)\right(3k+1)}{t}^{3k}]\\ \times [1+\frac{(k-1)(18k^4-91k^3+171k^2+187k+39)}{384(k+1{)}^3(2k+1\left)\right(3k+1\left)\right(4k+1)}{t}^{4k}]\cdots \end{array}$

For the expansion about the turning point, the nonintegral power in the argument of the hypergeometric function does not lead to singular expansion coefficients, but the resulting form is not conducive to useful results. Instead first apply an identity for Gauss hypergeometric functions,

$\begin{array}{l}{}_{2}F_1(a,b;c;z)=\frac{\Gamma \left(c\right)\Gamma (c-a-b)}{\Gamma (c-a)\Gamma (c-b)}{}_{2}F_1(a,b;a+b-c+1;1-z)\\ +\frac{\Gamma \left(c\right)\Gamma (a+b-c)}{\Gamma \left(a\right)\Gamma \left(b\right)}(1-z{)}^{c-a-b}{}_{2}F_1(c-a,c-b;c-a-b+1;1-z)\end{array}$

which in the present case becomes

$\begin{array}{l}{}_{2}F_1[\frac{1}{2},\frac{1}{k};\frac{1}{k}+1;x^k]=\frac{\sqrt{\pi }\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}{}_{2}F_1[\frac{1}{2},\frac{1}{k};\frac{1}{2};1-x^k]\\ -\frac{2}{k}\sqrt{1-x^k}{}_{2}F_1[\frac{1}{k}+\frac{1}{2},1;\frac{3}{2};1-x^k]\\ =\frac{\sqrt{\pi }\Gamma (\frac{1}{k}+1)}{\Gamma (\frac{1}{k}+\frac{1}{2})}\frac{1}{x}-\frac{2}{k}\sqrt{1-x^k}{}_{2}F_1[\frac{1}{k}+\frac{1}{2},1;\frac{3}{2};1-x^k]\end{array}$

Multiplying by a single power of x and recognizing the constant as the quarter period of the motion, the equation determining expansion coefficients is

$u=-\frac{2}{k}x\sqrt{1-x^k}{}_{2}F_1[\frac{1}{k}+\frac{1}{2},1;\frac{3}{2};1-x^k]$

The quantity $1-x^k$ is small and the hypergeometric function can be easily expanded with respect to that entire argument. Mathematica can handle expansion of the square root since all powers are integral there, so the evaluation proceeds in the same manner as the previous case. The result for an expansion depth of three terms beyond the initial unity is

$x_{\mathrm{turning}}\left(u\right)=[1-\frac{k}{4}{u}^2][1+\frac{(k-1)k^2}{96}{u}^4][1-\frac{(k-1)k^3(2k-11)}{2880}{u}^6]\cdots$

which when expanded to the same order is

$x_{\mathrm{turning}}\left(u\right)=1-\frac{k}{4}{u}^2+\frac{(k-1)k^2}{96}{u}^4-\frac{(k-1)k^3(4k-7)}{5760}{u}^6+\cdots$

These expansion coefficients again match the limits of multiple derivatives of powers of the rational function given above after being multiplied by the factorial of the power of the temporal variable in each term.

The expansions about the two endpoints should meet in the middle of the range with a judicious choice of number of terms. For the approximations given, the expansion about the turning point must be truncated by one term to allow meeting. An approximate inverse function can then be defined piecewise by reflection about the quarter period of the motion and the x-axis as necessary. It turns out to look like this:

This is a visualization of the motion defined by an inverse generalized hyperelliptic integral as a continuous analytic function of the exponent in the denominator of the integrand. Quite cool.

For a negative value of the exponent in the generalized hyperelliptic integral, first define a positive value p = −k. To avoid introducing an imaginary unit, the integral to be inverted will be taken to be

$\begin{array}{l}I(x,p)=\int \frac{dx}{\sqrt{\frac{1}{x^p}-1}}=\int dx\frac{x^{p/2}}{\sqrt{1-x^p}}\\ \phantom{I(x,p)}=\frac{2}{p+2}{x}^{(p+2)/2}{}_{2}F_1[\frac{1}{2},\frac{1}{p}+\frac{1}{2};\frac{1}{p}+\frac{3}{2};x^p]\end{array}$

where the integral is again evaluated with a Gauss hypergeometric function. Note that the second parameter of the Gauss hypergeometric function is $\frac{p+2}{2p}$ , a combination that links pairs of power potentials in both classical and quantum mechanics.

The quarter period for this inversion is

$\frac{T}{4}=\frac{2}{p+2}\frac{\Gamma (\frac{1}{p}+\frac{3}{2})\Gamma \left(\frac{1}{2}\right)}{\Gamma (\frac{1}{p}+1)\Gamma \left(1\right)}=\frac{\sqrt{\pi }\Gamma (\frac{1}{p}+\frac{1}{2})}{\Gamma \left(\frac{1}{p}\right)}$

so that the full period can be written equivalently

$T=4\sqrt{\pi }\frac{\Gamma (\frac{1}{p}+\frac{1}{2})}{\Gamma \left(\frac{1}{p}\right)}=2^{3-2/p}\pi \frac{\Gamma \left(\frac{2}{p}\right)}{\Gamma (\frac{1}{p}{)}^2}$

The full period for p = 1 is 2π because the Gauss hypergeometric function for that case is equivalent to an inverse sine plus the square root of a quadratic polynomial.

Lagrange inversion cannot here be used to indicate the general form of expansions due to the nonintegral power of the multiplicative factor. In order to determine the form of the expansion about the origin, set $x_{\mathrm{origin}}\left(t\right)=t^{2/(p+2)}f\left(t\right)$ in the equation to be inverted and again cancel a single power of the temporal variable:

$\frac{2}{p+2}[f\left(t\right){]}^{(p+2)/2}{}_{2}F_1[\frac{1}{2},\frac{1}{p}+\frac{1}{2};\frac{1}{p}+\frac{3}{2};t^{2p/(p+2)}\left[f\right(t){]}^p]=1$

By inspection $f\left(t\right)$ can only be a function of powers of $t^{2p/(p+2)}$ , so that the expansion about the origin is of the form

$x_{\mathrm{origin}}\left(t\right)=t^{2/(p+2)}\sum _{m=0}^{\infty }{c}_m{t}^{2mp/(p+2)}$

As a side note it now becomes clear that Lagrange inversion is not really needed to determine the form of expansions, but is worth including to see how well the method works in application. Lagrange inversion also provides a check on the accuracy of expansions for a positive exponent that will not be available for the negative exponent.

The expansion about the turning point is again a generalization of the behavior of sine about $\frac{\pi }{2}$ ,

$x_{\mathrm{turning}}\left(t\right)=1+\sum _{m=1}^{\infty }{d}_m[t-\frac{\sqrt{\pi }\Gamma (\frac{1}{p}+\frac{1}{2})}{\Gamma \left(\frac{1}{p}\right)}{]}^{2m}$

but with the appropriate quarter period. As for expansions for a positive exponent, it is more convenient to work with expansions in product form:

$\begin{array}{ll}{x}_{\mathrm{origin}}\left(t\right)=s^{2/(p+2)}\prod _{m=1}^{\infty }[1+a_m{s}^{2mp/(p+2)}]& s=\frac{p+2}{2}t\\ x_{\mathrm{turning}}\left(t\right)=\prod _{m=1}^{\infty }(1+b_m{u}^{2m})& u=t-\frac{\sqrt{\pi }\Gamma (\frac{1}{p}+\frac{1}{2})}{\Gamma \left(\frac{1}{p}\right)}\end{array}$

For the expansion about the origin, the method for determining the coefficient proceeds much as that for a positive exponent. The result for an expansion depth of four terms beyond the initial power is

$\begin{array}{l}{x}_{\mathrm{origin}}\left(s\right)=s^{2/(p+2)}[1-\frac{1}{(3p+2)}{s}^{2p/(p+2)}][1-\frac{(p+1)(p+2)}{2(3p+2{)}^2(5p+2)}{s}^{4p/(p+2)}]\\ \times [1-\frac{(p+1)(p+2)(2p^2+17p+6)}{3(3p+2{)}^3(5p+2\left)\right(7p+2)}{s}^{6p/(p+2)}]\\ \times [1-\frac{(p+1)(p+2)(162p^5+1881p^4+9323p^3+8494p^2+2804p+312)}{24(3p+2{)}^4(5p+2{)}^2(7p+2)(9p+2)}{s}^{8p/(p+2)}]\\ \times \cdots \end{array}$

For the expansion about the turning point, again apply an identity for Gauss hypergeometric functions,

$\begin{array}{l}{}_{2}F_1[\frac{1}{2},\frac{1}{p}+\frac{1}{2};\frac{1}{p}+\frac{3}{2};x^p]=\frac{\sqrt{\pi }\Gamma (\frac{1}{p}+\frac{3}{2})}{\Gamma (\frac{1}{p}+1)}{}_{2}F_1[\frac{1}{2},\frac{1}{p}+\frac{1}{2};\frac{1}{2};1-x^p]\\ -\frac{p+2}{p}\sqrt{1-x^p}{}_{2}F_1[\frac{1}{p}+1,1;\frac{3}{2};1-x^p]\\ =\frac{\sqrt{\pi }\Gamma (\frac{1}{p}+\frac{3}{2})}{\Gamma (\frac{1}{p}+1)}{x}^{-(p+2)/2}-\frac{p+2}{p}\sqrt{1-x^p}{}_{2}F_1[\frac{1}{p}+1,1;\frac{3}{2};1-x^p]\end{array}$

Multiplying by the remaining factors and simplifying, the equation determining expansion coefficients is

$u=-\frac{2}{p}{x}^{(p+2)/2}\sqrt{1-x^p}{}_{2}F_1[\frac{1}{k}+1,1;\frac{3}{2};1-x^p]$

The evaluation proceeds now as for the case of positive exponent values. The result for an expansion depth of three terms beyond the initial unity is

$x_{\mathrm{turning}}\left(u\right)=[1-\frac{p}{4}{u}^2][1+\frac{(p-3)p^2}{32}{u}^4][1-\frac{p^3(p^2-9p+17)}{192}{u}^6]\cdots$

which when expanded to the same order is

$x_{\mathrm{turning}}\left(u\right)=1-\frac{p}{4}{u}^2+\frac{(p-3)p^2}{32}{u}^4-\frac{(p-5)p^3(2p-5)}{384}{u}^6+\cdots$

An approximate inverse function is defined piecewise as for positive exponent values by matching the two expansions. The result looks like this:

For the Kepler potential p = 1 and nearby exponent values the expansion about the origin is quite good over the entire quarter period. The leading fractional power for negative exponent values allows one to construct a periodic motion even though the expansion itself does not have alternating signs.

The difference in behavior for positive and negative values of the exponent is most noticeable at the origin and half periods, resulting in two piecewise functions of distinctly different appearance. Once again, the complete code is available for extending the expansions by additional terms.

Uploaded 2015.05.22 — Updated 2017.07.31 analyticphysics.com